;;;; ----------------- monte carlo method -------------------

;;;---- monte-carlo method depends on two steps:
;;1. do many independent-repeated trials or samples(a trial or sample depends on a random-variable)
;;2. estimate (mean-value and error) samples been done

;;#!r6rs

;; In stream form, both trials and estimates are formed as streams
(library-directories '("." "/home/cl/gits/libs-for-chezscheme"))
(import (stream))

                                        ;aux: integers = (1 2 3 ...)
(define integers
  (let integers-from ([n 1])
    (stream-cons n (integers-from (+ n 1)))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;; monte-carlo
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; stream of samples
(define (sample-stream sampler); (sampler) => a sample of a random variable
  (stream-cons (sampler) (sample-stream sampler)))

;; sample-mean
(define (mean strm)
  (let ([sum-strm (stream-fold (lambda (acc x) (+ acc x))
                               0.0 strm)])
    (stream-map (lambda (n ele) (/ ele n))
                integers sum-strm)))

;; sample-variance
(define (variance strm)
  (let ([mean-strm (mean strm)]
        [mean-square-strm (mean (stream-map (lambda (x) (expt x 2)) strm))])
    (stream-map (lambda (n s mu)
                  (* (/ n (- n 1)) (- s (expt mu 2))))
                (stream-cdr integers)
                mean-square-strm
                mean-strm)))

;; stream of (mean-value . error), error can be: relative, absolute, sigma or kappa
(define (monte-carlo H sampler error-type)
  (let* ([samples (stream-map H (sample-stream sampler))]
         [mean-strm (mean samples)]
         [variance-strm (variance samples)]
         [method (lambda (error-type)
                   (lambda (n mu var)
                     (cons mu
                           (cond
                            [(eq? error-type 'relative) (* 1.96 (/ (sqrt (/ var n)) mu))]
                            [(eq? error-type 'absolute) (* 1.96 (sqrt (/ var n)))]
                            [(eq? error-type 'sigma) (sqrt var)]
                            (else (/ (sqrt (/ var n)) mu))))))]);;for kappa
    (stream-map (method error-type)
                integers
                mean-strm
                variance-strm))
  )


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;; Application

;;;;;;;;;;;;;;;;;;;;;;;
;; integration
;;;;;;;;;;;;;;;;;;;;;;;
;; problem: int( h(x), x in [0,1]^d )
(define (integral-with-errors func d)
  (let* ([sampler (lambda ()
                    (vector-map random (make-vector d 1.0)))])
    (monte-carlo (lambda (v)
                   (apply func (vector->list v)))
                 sampler
                 'absolute))
  )



#| ;;--examples--
;; (a) h = 2*x
;; precise result = 1
(let ([results (integral-with-errors (lambda (x) (* 2 x)) 1)])
  (stream-ref results 10000)
  )

;; (b) estimate Pi
;; use func H(x,y) = 4 (if x^2+y^2<1), or = 0 (others)
(let* ([H (lambda (x y) (if (< (+ (expt x 2) (expt y 2))
                               1)
                            4.0 0.0))]
       [Pi-stream (integral-with-errors H 2)]
       [N 10000])
  (stream-ref Pi-stream N)
  )
|#




